library(foreign)
## ---------------------------------------------------------------------------------------- ##
## Author: John Fox                                                                         ##
## Source: http://r.789695.n4.nabble.com/R-extend-summary-lm-for-hccm-td815004.html         ##
## Adapted by Tony Cookson.                                                                 ##
##        -- Only Change Made: Changed the name of the function (unwisely maybe)            ##
##           to summaryR from summaryHCCM.lm.  I also changed the spelling of consistent    ##
## ---------------------------------------------------------------------------------------- ##

summaryR.lm <- function(model, type=c("hc3", "hc0", "hc1", "hc2", "hc4"), ...){
  
  if (!require(car)) stop("Required car package is missing.")
  
  type <- match.arg(type)
  V <- hccm(model, type=type)
  sumry <- summary(model)
  table <- coef(sumry)
  table[,2] <- sqrt(diag(V))
  table[,3] <- table[,1]/table[,2]
  table[,4] <- 2*pt(abs(table[,3]), df.residual(model), lower.tail=FALSE)
  
  sumry$coefficients <- table
  p <- nrow(table)
  hyp <- cbind(0, diag(p - 1))
  sumry$fstatistic[1] <- linearHypothesis(model, hyp,white.adjust=type)[2,"F"]
  
  print(sumry)
  cat("Note: Heteroscedasticity-consistent standard errors using adjustment", type, "\n")
  
}


dat<-read.dta(file="POQdataset-ts-additions.dta")

dat$turnover<-as.numeric(dat$year%in%c(1953, 1961, 1963, 1969, 1974, 1977, 1981, 1989, 1993, 2001))
attach(dat)


# basic plots
dat1<-subset(dat, lapp2<50)
plot(app2~unempl, data=dat1, ylim=c(20, 90), xlim=c(2, 12))
abline(lm(app2~unempl, data=dat1))
dat2<-subset(dat, lapp2>=50)
points(app2~unempl, data=dat2, col="gray40")
abline(lm(app2~unempl, data=dat2), col="gray40")

model.bg<-lm(app2~lapp2+growth+unempl+d_Infl_new+logKcas+logVcas+logAIcas)
summaryR.lm(model.bg)

lapp2Xgrowth=lapp2*growth
model.sd.1<-lm(app2~lapp2+growth+lapp2Xgrowth+unempl+d_Infl_new+logKcas+logVcas+logAIcas)
summaryR.lm(model.sd.1)

lapp2Xunempl=lapp2*unempl
model.sd.2<-lm(app2~lapp2+growth+unempl+lapp2Xunempl+d_Infl_new+logKcas+logVcas+logAIcas)
summaryR.lm(model.sd.2)

model.sd.3<-lm(app2~lapp2+growth+lapp2Xgrowth+unempl+lapp2Xunempl+d_Infl_new+logKcas+logVcas+logAIcas)
summaryR.lm(model.sd.3)

AIC(model.bg, model.sd.1, model.sd.2, model.sd.3)
BIC(model.bg, model.sd.1, model.sd.2, model.sd.3)

model.sd.2a<-lm(app2~lapp2+growth+unempl+lapp2Xunempl+d_Infl_new+logKcas+logVcas+logAIcas+turnover)
summaryR.lm(model.sd.2a)


# drawing the coefficients out of the asymptotic VCV
library(mvtnorm)
beta.draws<-rmvnorm(1000, mean=coefficients(model.sd.2), sigma=hccm(model.sd.2))


# collect and plot the MEs for each value of y_lag
y.seq<-seq(from=23, to=87, by=1)
me<-matrix(data=NA, ncol=3, nrow=length(y.seq))

for(i in 1:length(y.seq)){
  
  y.test<-y.seq[i]
  me.draws<-beta.draws[,4]+beta.draws[,5]*y.test
  me[i,]<-quantile(me.draws, probs=c(0.05, 0.5, 0.95))
  
  
}

par(mar=c(5, 4.5, 4, 2))
plot(me[,2]~y.seq, type="l", ylim=c(-4, 2), xlim=c(20,90), main=c("Instantaneous Marginal Effect", "of Unemployment on Pres. Approval"), xlab=expression(Approval[(t-1)]), ylab=expression(partialdiff*Approval[t]/partialdiff*Unemployment[t]))
lines(me[,1]~y.seq, lty=2)
lines(me[,3]~y.seq, lty=2)
abline(h=0, lty=3)
legend("bottomleft", lty=c(1,2), c("marginal effect", "90% confidence interval"))

dev.copy2eps(file="geys_instant_mfx_r.eps")




unempl.seq<-seq(from=2.5, to=10.5, by=0.1)
me<-matrix(data=NA, ncol=3, nrow=length(unempl.seq))

for(i in 1:length(unempl.seq)){
  
  unempl.test<-unempl.seq[i]
  me.draws<-beta.draws[,2]+beta.draws[,5]*unempl.test
  me[i,]<-quantile(me.draws, probs=c(0.05, 0.5, 0.95))
  
  
}

par(mar=c(5, 4.5, 4, 2))
plot(me[,2]~unempl.seq, type="l", xlim=c(2,11), ylim=c(0.4, 1.2), main=c("Instantaneous Marginal Effect", "of Lag on Presidential Approval"), xlab=expression(Unemployment[t]), ylab=expression(partialdiff*Approval[t]/partialdiff*Approval[(t-1)]))
lines(me[,1]~unempl.seq, lty=2)
lines(me[,3]~unempl.seq, lty=2)
abline(h=0, lty=3)
legend("bottomleft", lty=c(1,2), c("marginal effect", "90% confidence interval"))

dev.copy2eps(file="geys_instant_mfx_lag.eps")





# plot long term trajectories for the Geys model

x.8<-matrix(dat=c(1,50,3.366683,8,8*50,-.0182841,.5229144,1.25769,.5709695), nrow=1)
y.8<-rbind(50,x.8%*%t(beta.draws))

for(i in 3:9){
  
  y.8<-rbind(y.8,NA)
  for(j in 1:1000){
    x.8<-matrix(dat=c(1,y.8[i-1,j],3.366683,8,8*y.8[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.8[i,j]<-x.8%*%t(t(beta.draws[j,]))
  }
  
}


x.7<-matrix(dat=c(1,50,3.366683,7,7*50,-.0182841,.5229144,1.25769,.5709695), nrow=1)
y.7<-rbind(50,x.7%*%t(beta.draws))

for(i in 3:9){
  
  y.7<-rbind(y.7,NA)
  for(j in 1:1000){
    x.7<-matrix(dat=c(1,y.7[i-1,j],3.366683,7,7*y.7[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.7[i,j]<-x.7%*%t(t(beta.draws[j,]))
  }
  
}

x.4<-matrix(dat=c(1,50,3.366683,4,4*50,-.0182841,.5229144,1.25769,.5709695), nrow=1)
y.4<-rbind(50,x.4%*%t(beta.draws))

for(i in 3:9){
  
  y.4<-rbind(y.4,NA)
  for(j in 1:1000){
    x.4<-matrix(dat=c(1,y.4[i-1,j],3.366683,4,4*y.4[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.4[i,j]<-x.4%*%t(t(beta.draws[j,]))
  }
              
}

x.3<-matrix(dat=c(1,50,3.366683,3,3*50,-.0182841,.5229144,1.25769,.5709695), nrow=1)
y.3<-rbind(50,x.3%*%t(beta.draws))

for(i in 3:9){
  
  y.3<-rbind(y.3,NA)
  for(j in 1:1000){
    x.3<-matrix(dat=c(1,y.3[i-1,j],3.366683,3,3*y.3[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.3[i,j]<-x.3%*%t(t(beta.draws[j,]))
  }
  
}

y.ch=y.4-y.3

t<-1:8
traj.quant<-apply(X=y.ch[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
y.3.traj<-apply(X=y.3[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
y.4.traj<-apply(X=y.4[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))

plot(traj.quant[2,]~t, type="l", xlab="time")
lines(traj.quant[1,]~t, lty=2)
lines(traj.quant[3,]~t, lty=2)
t.alt<-t+0.25
library(gplots)
plotCI(x=t, y=traj.quant[2,], ui=traj.quant[3,], li=traj.quant[1,], ylab="difference in approval rating", xlab="time (in quarters)", main = c("simulated difference in approval trajectory", "with 90% CI error bars"), gap=0.5, xlim=c(0,8), ylim=c(-12,0.1), xaxp=c(0,8,8))
points(0,0)
dev.copy2eps(file="approval-dif-trajectory.eps")

plotCI(x=t, y=y.3.traj[2,], ui=y.3.traj[3,], li=y.3.traj[1,], ylab="approval rating", xlab="time (in quarters)", main = c("simulated approval trajectory", "with 90% CI error bars"), gap=0.5, xlim=c(0,8.22), ylim=c(45, 76), xaxp=c(0,8,8))
points(0, 50)
t.alt<-t+0.25
plotCI(x=t.alt, y=y.4.traj[2,], ui=y.4.traj[3,], li=y.4.traj[1,], gap=0.5, col="gray45", add=T)
points(0.25, 50, col="gray45")
legend("topleft", lty=c(1,1), col=c("black", "gray45"), legend=c("3% unemployment", "4% unemployment"))
dev.copy2eps(file="approval-trajectory.eps")


y.7.traj<-apply(X=y.7[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
plotCI(x=t, y=y.7.traj[2,], ui=y.7.traj[3,], li=y.7.traj[1,], ylab="approval rating", xlab="time (in quarters)", main = c("simulated approval trajectory", "with 90% CI error bars"), gap=0.5, xlim=c(0, 8.22), ylim=c(32, 52), xaxp=c(0,8,8))

y.8.traj<-apply(X=y.8[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
plotCI(x=t.alt, y=y.8.traj[2,], ui=y.8.traj[3,], li=y.8.traj[1,], ylab="approval rating", xlab="time (in quarters)", main = c("simulated approval trajectory", "with 90% CI error bars"), gap=0.5, col="gray45", add=T)

points(0, 50)
points(0.25, 50, col="gray45")

legend("bottomleft", lty=c(1,1), col=c("black", "gray45"), legend=c("7% unemployment", "8% unemployment"))
dev.copy2eps(file="approval-trajectory-highU.eps")


y.ch=y.8-y.7
traj.quant<-apply(X=y.ch[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
library(gplots)
plotCI(x=t, y=traj.quant[2,], ui=traj.quant[3,], li=traj.quant[1,], ylab="difference in approval rating", xlab="time (in quarters)", main = c("simulated difference in approval trajectory", "with 90% CI error bars"), gap=0.5, xlim=c(0,8), ylim=c(-5, 0.1), xaxp=c(0,8,8))
points(0,0)
dev.copy2eps(file="approval-dif-trajectory-highU.eps")




# simulate 8% unemployment for 4 quarters, then 6% unemployment for 4 quarters

x.int<-matrix(dat=c(1,50,3.366683,8,8*50,-.0182841,.5229144,1.25769,.5709695), nrow=1)
y.int<-rbind(50,x.int%*%t(beta.draws))

for(i in 3:5){
  
  y.int<-rbind(y.int,NA)
  for(j in 1:1000){
    x.int<-matrix(dat=c(1,y.int[i-1,j],3.366683,8,8*y.int[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.int[i,j]<-x.int%*%t(t(beta.draws[j,]))
  }
  
}

for(i in 6:9){
  
  y.int<-rbind(y.int,NA)
  for(j in 1:1000){
    x.int<-matrix(dat=c(1,y.int[i-1,j],3.366683,6,6*y.int[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.int[i,j]<-x.int%*%t(t(beta.draws[j,]))
  }
  
}

y.int.traj<-apply(X=y.int[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
plotCI(x=t, y=y.int.traj[2,], ui=y.int.traj[3,], li=y.int.traj[1,], ylab="approval rating", xlab="time (in quarters)", main = c("simulated approval trajectory", "with 90% CI error bars"), gap=0.5, xlim=c(0, 8.22), ylim=c(32, 56), xaxp=c(0,8,8))
points(0,50)
abline(v=4.5, lty=2)
text(x=2, y=35, labels="8% unemployment")
text(x=6.5, y=34, labels="6% unemployment")


for(i in 6:9){
  
  y.int<-rbind(y.int,NA)
  for(j in 1:1000){
    x.int<-matrix(dat=c(1,y.int[i-1,j],3.366683,4,4*y.int[i-1,j],-.0182841,.5229144,1.25769,.5709695), nrow=1)
    y.int[i,j]<-x.int%*%t(t(beta.draws[j,]))
  }
  
}

t.alt<-t+0.25
y.int.traj<-apply(X=y.int[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
plotCI(x=t.alt[5:8], y=y.int.traj[2,5:8], ui=y.int.traj[3,5:8], li=y.int.traj[1,5:8], gap=0.5, add=T, col="gray45")
text(x=6.5, y=36, labels="4% unemployment", col="gray45")

dev.copy2eps(file="approval-improving-economy.eps")






# simulate trajectories with and without 1000 casualties at 4% and 8% unemployment
x.4.cas<-matrix(dat=c(1,50,3.366683,4,4*50,-.0182841,0,log(1000),0), nrow=1)
y.4.cas<-rbind(50,x.4.cas%*%t(beta.draws))

for(i in 3:9){
  
  y.4.cas<-rbind(y.4.cas,NA)
  for(j in 1:1000){
    x.4.cas<-matrix(dat=c(1,y.4.cas[i-1,j],3.366683,4,4*y.4.cas[i-1,j],-.0182841,0,log(1000),0), nrow=1)
    y.4.cas[i,j]<-x.4.cas%*%t(t(beta.draws[j,]))
  }
  
}

y.4.cas.traj<-apply(X=y.4.cas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))

x.4.nocas<-matrix(dat=c(1,50,3.366683,4,4*50,-.0182841,0,0,0), nrow=1)
y.4.nocas<-rbind(50,x.4.nocas%*%t(beta.draws))

for(i in 3:9){
  
  y.4.nocas<-rbind(y.4.nocas,NA)
  for(j in 1:1000){
    x.4.nocas<-matrix(dat=c(1,y.4.nocas[i-1,j],3.366683,4,4*y.4.nocas[i-1,j],-.0182841,0,0,0), nrow=1)
    y.4.nocas[i,j]<-x.4.nocas%*%t(t(beta.draws[j,]))
  }
  
}

y.4.nocas.traj<-apply(X=y.4.nocas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))

y.ch=y.4.cas-y.4.nocas
traj.quant<-apply(X=y.ch[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
library(gplots)
par(mar=c(5, 4, 4, 2))
plotCI(x=t, y=traj.quant[2,], ui=traj.quant[3,], li=traj.quant[1,], ylab="(approval|1000 casualties) - (approval|0 casualties)", xlab="time (in quarters)", main = c("simulated change in approval over time", "as Vietnam Casualties rise, with 90% CI error bars"), gap=0.5, xlim=c(0,8.22), ylim=c(-30, 0.1), xaxp=c(0,8,8))
points(0,0)

# simulate trajectories with and without 1000 casualties
x.8.cas<-matrix(dat=c(1,50,3.366683,8,8*50,-.0182841,0,log(1000),0), nrow=1)
y.8.cas<-rbind(50,x.8.cas%*%t(beta.draws))

for(i in 3:9){
  
  y.8.cas<-rbind(y.8.cas,NA)
  for(j in 1:1000){
    x.8.cas<-matrix(dat=c(1,y.8.cas[i-1,j],3.366683,8,8*y.8.cas[i-1,j],-.0182841,0,log(1000),0), nrow=1)
    y.8.cas[i,j]<-x.8.cas%*%t(t(beta.draws[j,]))
  }
  
}

y.8.cas.traj<-apply(X=y.8.cas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))

x.8.nocas<-matrix(dat=c(1,50,3.366683,8,8*50,-.0182841,0,0,0), nrow=1)
y.8.nocas<-rbind(50,x.8.nocas%*%t(beta.draws))

for(i in 3:9){
  
  y.8.nocas<-rbind(y.8.nocas,NA)
  for(j in 1:1000){
    x.8.nocas<-matrix(dat=c(1,y.8.nocas[i-1,j],3.366683,8,8*y.8.nocas[i-1,j],-.0182841,0,0,0), nrow=1)
    y.8.nocas[i,j]<-x.8.nocas%*%t(t(beta.draws[j,]))
  }
  
}

y.ch=y.8.cas-y.8.nocas
traj.quant<-apply(X=y.ch[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
library(gplots)
t.alt=t+0.25
plotCI(x=t.alt, y=traj.quant[2,], ui=traj.quant[3,], li=traj.quant[1,], gap=0.5, add=T, col="gray45")
points(0.25,0,col="gray45")
legend("bottomleft", lty=c(1,1), pch=c(1,1), col=c("black","gray45"), legend=c("4% unemployment", "8% unemployment"))
dev.copy2eps(file="approval-dif-cas-trajectory.eps")



traj.8.cas<-apply(X=y.8.cas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
plotCI(x=t, y=traj.8.cas[2,], ui=traj.8.cas[3,], li=traj.8.cas[1,], gap=0.5, add=F, col="gray45", xlim=c(0,8.22), ylim=c(20, 70), xaxp=c(0,8,8), main="8% unemployment", xlab="time (in quarters)", ylab="approval rating")
points(0,50)

traj.8.nocas<-apply(X=y.8.nocas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
t.alt=t+0.25
plotCI(x=t.alt, y=traj.8.nocas[2,], ui=traj.8.nocas[3,], li=traj.8.nocas[1,], gap=0.5, add=T, col="black")
points(0.25,50,col="gray45")

legend("topleft", lty=c(1,1), pch=c(1,1), col=c("black","gray45"), legend=c("no casualties", "1000 casualties"))

dev.copy2eps(file="approval-8pct-trajectory.eps")



traj.4.cas<-apply(X=y.4.cas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
t.alt=t+0.25
plotCI(x=t.alt, y=traj.4.cas[2,], ui=traj.4.cas[3,], li=traj.4.cas[1,], gap=0.5, add=F, col="gray45", xlim=c(0,8.22), ylim=c(20, 70), xaxp=c(0,8,8), main="4% unemployment", xlab="time (in quarters)", ylab="approval rating")
points(0,50)

traj.4.nocas<-apply(X=y.4.nocas[2:9,], MARGIN=1, FUN=quantile, probs=c(0.05, 0.5, 0.95))
t.alt=t+0.25
plotCI(x=t.alt, y=traj.4.nocas[2,], ui=traj.4.nocas[3,], li=traj.4.nocas[1,], gap=0.5, add=T, col="black")
points(0.25,50,col="gray45")

legend("topleft", lty=c(1,1), pch=c(1,1), col=c("black","gray45"), legend=c("no casualties", "1000 casualties"))


dev.copy2eps(file="approval-4pct-trajectory.eps")